Assinaturas hidrológicas para bacias hidrográficas do SIN

Pré-requisitos

Pacotes necessários:

#if(!require(PSF)) install.packages("PSF")
#if(!require(timetk)) remotes::install_github("business-science/timetk")
if(!require(EflowStats)) remotes::install_github("USGS-R/EflowStats")


pacotes <- c(
  "here",
  "usethis",
  "data.table",
  "HEobs",
#  "PSF",
  "tidyverse",
  "lubridate",
  "fs",
  "checkmate",
#  "xts",
#  "hydroGOF",
#  "ModelMetrics",
#  "forecast",
#  "timetk",
  "EflowStats",
  "hydrostats",
 #"NbClust",
 #  "cluster",  
 #  "cluster.datasets", 
  "cowplot", 
  # "clValid",
  "ggfortify", 
  #"clustree",
  #"dendextend",
  #"factoextra",
  #"FactoMineR",
  #"corrplot",
  #"GGally",
  #"ggiraphExtra",
  "kableExtra",
  "tidymodels",
  "fasstr"
)
# Carregar os pacotes
easypackages::libraries(pacotes)

Scripts:

source(here('R', 'load-data.R'))
source(here('R', 'utils.R'))

Dados de vazão

Os dados importados de vazão devem ser regularmente espaçados no tempo. Esta adequação das séries diárias, se necessária, pode ser realizada com a função complete_dates() do pacote {lhmetools}. Assim assegura-se que todas estações possuam 1 observação por dia e sem datas faltantes.

Metadados

data_link <- "https://www.dropbox.com/s/d40adhw66uwueet/VazoesNaturaisONS_D_87UHEsDirceuAssis_2018.dat?dl=1"
qnat_meta <- extract_metadata(NA_character_, informative = TRUE)
glimpse(qnat_meta)

Dados

qnat_data <- qnat_dly_ons() %>%
  select(date, qnat, code_stn) %>%
  lhmetools::complete_dates(group = "code_stn")
glimpse(qnat_data)

# Incluindo nome da estacao
qnat_data <- qnat_data %>%
  full_join(select(qnat_meta, estacao_codigo, nome_estacao),
            by = c(code_stn = "estacao_codigo")
            ) %>%
  rename("name_stn" = nome_estacao)

Assinaturas hidrológicas para 1 posto

Preparação dos dados para aplicação da função calc_magnifSeven com a opção de ano hidrológico (water year) para o posto 74 da ONS (G. B. Munhoz).

qnat_posto <- qnat_data %>% 
  sel_station(.,station = 74)  
glimpse(qnat_posto)
summary(qnat_posto)

magnif7(stn_data = select(qnat_posto, date, qnat))

Dados agrupados por postos

by_stn <- qnat_data %>% 
  group_by(code_stn) %>%
  nest()

Gráfico da sazonalidade da fração de vazão anual.

Fração mensal da vazão anual

# https://jcoliver.github.io/learn-r/008-ggplot-dendrograms-and-heatmaps.html
# glimpse(by_stn[["data"]][[1]])

q_mon_clim <- by_stn %>%
  unnest(cols = "data") %>%
  group_by(code_stn, month = lubridate::month(date)) %>%
  select(-name_stn) %>%
  summarise(q_med = mean(qnat, na.rm = TRUE))

q_ann_clim <- q_mon_clim %>%
  summarise(q_tot = sum(q_med))

q_mon_clim <- q_mon_clim %>% 
  full_join(q_ann_clim) %>%
  mutate(q_frac = q_med/q_tot * 100) %>%
  ungroup() %>%  
  mutate(code_stn = as.factor(code_stn),
         month = as.factor(month))
  
psych::describe(q_mon_clim) %>% 
  relocate(skew, kurtosis)

Agrupamento dos postos

qlong <- select(q_mon_clim, code_stn, month, q_frac) %>%
  pivot_wider(names_from = "month", 
              values_from = "q_frac",
              names_prefix = "qfrac_"
              )
qlong_scaled <- qlong 
qlong_scaled[,-1] <- scale(qlong_scaled[,-1])

# Run clustering
qmat <- as.matrix(qlong_scaled[, -1])
rownames(qmat) <- qlong_scaled$code_stn

clustd <- dist(x = qmat)
cd <- hclust(d = clustd, method="ward.D2")
#cd <- hclust(d = dist(x = qmat))
# optclus <- sapply(2:20, 
#                    function(i) 
#                      summary(cluster::silhouette(cutree(cd, k = i), clustd))$avg.width
#                    )
# optclus
# which.max(optclus) # 2
q_dendro <- as.dendrogram(cd)

# Create dendro
dendro_plot <- ggdendro::ggdendrogram(data = q_dendro, rotate = TRUE)+
   theme(axis.text.y = element_text(size = 8))

# Preview the plot
dendro_plot

# linhas da ordem dos postos
q_order <- order.dendrogram(q_dendro)
# ordem dos postos
as.integer(as.vector(qlong_scaled$code_stn[q_order]))

Heat map da fração mensal da vazão anual com postos ordenados pelo agrupamento.

q_mon_clim_trans <-  q_mon_clim %>%
  mutate(code_stn = factor(x = code_stn,
                               levels = qlong_scaled$code_stn[q_order], 
                               ordered = TRUE),
         q_frac = round((q_frac/2))*2
         )

# cbind(t = q_mon_clim_trans$q_frac, o = q_mon_clim$q_frac)
cols <- c("firebrick1", "lightpink3", "lightsteelblue3", 
          "cadetblue1", "cornflowerblue","blue", "navyblue")

q_mon_clim_trans %>%
  
  ggplot(aes(x = month, y = code_stn)) +
  geom_tile(aes(fill = q_frac), colour = "white") +
  #geom_tile(aes(fill = factor(q_frac)), colour = "white") +
  scale_x_discrete(expand = c(0,0)) +
  theme_bw() +
  #scale_fill_viridis_c()
  #scale_fill_viridis_c(guide = "legend")
  scale_fill_gradientn( "Vazão (% da média anual)",
                       colours = cols,
                       #guide = "legend",
                       breaks = seq(0, 20, by = 2)
                       ) +
  #scale_fill_manual(values = cols) +
  theme(legend.title = element_text(angle = 90, vjust = 1),
        legend.key.height = unit(1.5, units = "cm")
        )

    # scale_fill_distiller(
    #                      palette='RdYlBu', 
    #                      direction = 1,
    #                      breaks = seq(0, 20, by = 2),
    #                      #guide = "legend"
    #                     )
  

library(heatmaply)
index_cols <- c(8:12, 1:7) + 1
qdata <- as.data.frame(qlong[, index_cols])
names(qdata) <- tolower(month.abb[index_cols-1])

lut_postos <- qnat_data %>% 
  select(contains("stn")) %>% 
  distinct() 
postos <- lut_postos$name_stn
postos <- setNames(postos, lut_postos$code_stn)
  
rownames(qdata) <- postos[as.character(qlong$code_stn)]

heatmaply(qdata, 
          k_row = 4, 
          Colv = FALSE,
          #k_col = 3, 
          #scale = "column"
          scale_fill_gradient_fun = scale_fill_gradientn( "Vazão \n (% da média anual)",
                       colours = cols,
                       #guide = "legend",
                       breaks = seq(0, 20, by = 2)
                       ),
          fontsize_row = 5
          #seriate = "mean",
          #seriate = "OLO"
          #seriate = "GW"
          #seriate = "none"
          )

Testando fasstr

https://bcgov.github.io/fasstr/articles/fasstr.html

Será relevante considerar a diferenças nos períodos dos anos hidrológicos por posto?

check_season_plots <- qnat_data %>%
  rename("Date" = date) %>%
  fasstr::plot_daily_stats(
    values = "qnat",
    groups = "code_stn",
    start_year = 1991,
    end_year = 2010,
    log_discharge = TRUE,
    add_year = 2001, 
    complete_years = TRUE, 
    include_title = TRUE
  )

codes <- names(check_season_plots) %>%
  readr::parse_number()

lookup <- qnat_data %>% 
  select(contains("stn")) %>% 
  distinct() %>%
  slice(order(codes)) %>%
  arrange(code_stn)

plot_l <- map(seq_along(codes), 
    function(i) {
      # i = 1
      nm <- filter(lookup, code_stn == codes[i]) %>%
        pull(name_stn)
      check_season_plots[[i]] + ggtitle(paste0("Posto: ", nm))
    })
plot_l

Assinaturas hidrológicas para todos postos

seven_stats <- by_stn %>%
  #ungroup() %>%
  mutate(stats = purrr::map(data, ~.x %>% select(date, qnat) %>% magnif7))

seven_stats_tidy <- seven_stats %>%
  select(stats) %>%
  unnest(cols = "stats") %>%
  pivot_wider(names_from = indice, values_from = statistic) %>%
  ungroup()
seven_stats_tidy
saveRDS(seven_stats_tidy, file = here("output", "seven_stats.RDS"))